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Abstract 

A nonstandard application of bivariate polynomial interpolation is discussed: the 
implicitization of a rational algebraic curve given by its parametric equations. Three 
different approaches using the same interpolation space are considered, and their 
respective computational complexities are analyzed. Although the techniques em¬ 
ployed are usually associated to numerical analysis, in this case all the computations 
are carried out using exact rational arithmetic. The power of the Kronecker product 
of matrices in this application is stressed. 
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1 Introduction 


Curve implicitization, which consists of hnding the implicit equation of a curve 
C given by a rational parametrization, is an important problem in computer 
aided geometric design, and several theoretical results which help to its so¬ 
lution have been developed in the helds of classical algebraic geometry and 
computer algebra. Some of the methods for the effective computation of the 
implicit equation are based on interpolation. 

We will consider three types of interpolation problems which share the same 
interpolation space. The different ways of choosing the interpolation nodes will 
lead to linear systems with very different coefficient matrices: an unstructured 

^ Corresponding author. E-mail: ana.marco@uah.es 
^ E-mail: jjavier.martinez@uah.es 


Preprint submitted to Elsevier Preprint 


2 February 2008 





one, the transpose of a Vandermonde matrix, and the Kronecker product of 
two Vandermonde matrices. 

The Kronecker product structure of the matrix in the third method will make 
much less expensive the computational cost of the process [12, 22] and will 
introduce a high degree of parallelism. 

In this application of interpolation exact arithmetic is used, so stability is not 
as important as efficiency and we do not have to worry about matters which 
are very important in numerical analysis such as the ill-conditioning of Van¬ 
dermonde matrices [9]. In fact, the algorithms to be described use techniques 
from numerical linear algebra but all the computations are carried out in ex¬ 
act arithmetic, which allows them to be applied with no difficulty even when 
polynomials of very high degree are involved. 

Finally, the function to be interpolated is itself a bivariate polynomial, so 
there will be no interpolation error, and consequently in this case the use of 
interpolation will not be related to approximation theory. 


In order to make our exposition as clear and complete as possible, we will 
begin with a very small example which will be useful to illustrate the different 
approaches to the problem. Let us consider the hyperbola with parametric 
equations 


{x{t),y{t)) 


1 + t 3 + t\ 

2 + f’ 4 +J' 


Its implicit equation is 

2 — 3y — X + 2xy = 0. 


In the following sections it will be shown how this implicit equation can be 
obtained. 


The hrst problem to be addressed is the choice of an appropriate interpolation 
space. Let us observe that the implicit equation of any degree n polynomial 
rational parametric curve, with the same denominator in both components of 
the parametrization, is a degree n algebraic curve (see Section 15.4 of [20]). 
However, this is not necessarily true when two different denominators are 
considered in the parametrization. In our example, both parametric equations 
have degree 1 while the implicit equation has total degree 2. This is the reason 
why, when considering the general situation of two different denominators, it 
is more natural to consider the coordinate degree (i.e. the degree in x, the 
degree in y) instead of the total degree. The following result (see [15] and [21]) 
gives us the precise form of the implicit equation. 


Theorem 1. Let P = ix{t) = y^,|/(^) = ® proper rational 


parametrization of the irreducible curve C defined by F{x,y), and let 
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gcd{ui,vi) = gcd{u 2 ,V 2 ) = 1. Then max{degt{ui), degt{vi)} = degy{F) and 
max{degt{u 2 ),degt{v 2 )} = deg^{F). 

Theorem 1 tells us that the polynomial F{x,y) dehning the implicit equa¬ 
tion of the curve C belongs to the polynomial space where m = 

max{degt{u 2 ),degt{v 2 )} and n = max{degt{ui),degt{vi)}. The dimension of 
nm,n(a^, y) is = (m -|- l)(n -|- 1), and a basis is given by 

{xV = 0, • ■ ■ , m; j = 0, • ■ ■ , n} = 

= {1, y, ■ ■ ■ , y"", x,xy,--- , xy”, • • ■ , x”^, x^y, • ■ ■ , x”^?/”}. 

Moreover degx{F{x,y)) = m and degy{F{x,y)) = n, and therefore there is no 
interpolation space 11^,5 (x,?/) with r < m oi s < n such that F{x,y) belongs 
to n^, 5 (x,?/). 

In addition, the selection of the interpolation space Ylrn,n{,x, y) is also snitable 
because in practice the implicit representation of a rational parametric curve 
is a dense polynomial [13]. 

So in the application of interpolation we are considering the interpolation space 
is given and we have to choose appropriate Lagrange interpolation nodes. As 
we will see, for obtaining a linear system with a structnred coefficient matrix 
we will need the introduction of resultants, a tool which is widely used in 
computer algebra and has a variety of applications [4]. 

In the next three sections we will discuss three different approaches to the 
implicitization problem by using interpolation, and in Section 5 we will analyze 
their respective compntational complexities. 


2 An interpolation problem with an unstructured coefficient ma¬ 
trix 


The first approach, included in [11] along with other approaches to impliciti¬ 
zation, leads to a non-structured interpolation problem whose solution is not 
nnique. 

Since the dimension of the linear space of solutions will be seen to be al¬ 
most snrely 1, the implicit equation can be computed in the following way. 
As we know the interpolation space (it is Ilm,n{x,y)), we may evaluate the 
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parametrization P = {x{t),y{t)) at some finite set of values of t (for example 
t = 0,1, 2,...) for obtaining N distinct rational interpolation nodes 

Vi) = ixiti),yiti)) : i = 1,..., iV}. 


Then we can formulate an interpolation problem using those interpolation 
nodes, in this case with all the interpolation data equal to zero. So, we have 
an interpolation problem in which the corresponding linear system is a non- 
structured homogeneous one, and we know from Section 1 that a nontrivial 
solution always exists. 

If we denote by A the coefficient matrix of this homogeneous linear system 
and r = rank{A) then the set of all solutions of Ax = 0 is the nullspace of A, 
and its dimension N — r is great or equal than 1. 

Let us recall that Bezout’s theorem states that two plane eurves of degree m+n 
without common components have at most (m + n)^ common complex points. 
As a consequence of this theorem, the probability of having rank{A) < iV—1 is 
negligible, since we are prescribing N = (m+l)(n+l) interpolation conditions 
at rational points. However, if that happens we can add a new interpolation 
condition to hnd the correct implicit equation. 


In our example the rank is necessarily equal to 3 because m = n = 1 and 
therefore (m + n)^ = 4 = (m + l)(n + 1). If we consider t = 0,1,2,3 the 
interpolation nodes in the corresponding order are 


1 3 
2’4 


2 4 
3’5 


3 5 
4’6 


4 6 
5’7 


and the coefficient matrix of the homogeneous linear system Ac = 0 is 

^ 3 1 3 ^ 

■^428 


A = 


1 4 2 _8_ 

5 3 15 

1 5 3 5 
■^648 

1 6 4 ^ 

7 5 35y 


The solution vector is 


c=(2,-3,-l,2f 

or any multiple of it, and so the implicit equation of the given curve is 


2 — 3y — X + 2xy = 0. 
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3 The use of resultants and the Vandermonde matrix 


The computation of the implicit equation of a rationally parametrized curve 
can be carried out by computing the resultant of two polynomials. The fol¬ 
lowing theorem provides the way of doing it (see [15] and [21]). 


Theorem 2. Let P = (x{t) = j^,y{t) = be a proper rational 


parametrization of an irredueible eurve C, with gcd{ui,Vi) = gcd{u2,V2) = 1. 
Then the polynomial defining C is Resfiuiit) — xvi{t), U 2 {t) — yv 2 (t)) (the re¬ 
sultant with respect to t of the polynomials Ui(t) — xvi(t) and U 2 (t) — yv 2 {t)). 


Taking this theorem into account, the polynomial defining the implicit equa¬ 
tion of the curve of the example introduced in Section 1 is 

Resfiifi -ft) — x(2 -\- 1), (3 -|- t) — 1/(4 -1- t)), 


which is precisely 

2 — ?,y — X + 2xy. 

The resultant of two polynomials can be computed, for example, by using the 
command resultant of the symbolic computation system Maple. 

However, let us point out here that, in general, the computation of Resfiufit) — 
xviit), U 2 (t) —yv 2 (t)) is not a trivial task. It is the determinant of the Sylvester 
(or Bezout) matrix of ufit) — xvi{t) and U 2 (t) — yv 2 (t) (see [20], for example), 
and therefore its computation requires the expansion of the determinant of a 
matrix whose entries are polynomials in the variables x and y. The symbolic 
expansion of a determinant is a computer algebra problem which requires a lot 
of time and space to be solved due to the problem of intermediate expression 
swell. This fact is recognized in [13] where it is said that the bottleneck of the 
algorithm for implicitizing rational surfaces is the symbolic expansion of the 
determinant. As it can be read in [5], one of the most interesting approaches 
for the symbolic expansion of the determinant is based on interpolation, and 
it is the one presented in [13] and [14]. 

In that approach the N interpolation data are obtained by evaluating the 
Sylvester (or Bezout) matrix at the N interpolation nodes 

{{Pi,P 2 ) ■ k = 0,..., N — 1; Pi,P 2 distinct primes}, 

and computing the corresponding constant determinants. This clever selection 
of the nodes reduces the solution of the interpolation problem to the solution 
of a linear system of order N whose coefficient matrix is the transpose of a 
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nonsingular (since different monomials evaluate to different values) Vander¬ 
monde matrix (see, for example, [17] for the expression of the Vandermonde 
matrix). 


We illustrate this approach with our example. 


We start by showing the Sylvester matrix of the polynomials p{t) = 1 -|- t — 
x{2 + t) and q{t) = 3 -|- t — ?/(4 -|- t), which is 


S 


\ — X 1 — 2x 
- 1/ 3 - % 


In order to compute its determinant by using this interpolation approach, we 
consider pi = 2 and p 2 = 3, and so the iV = 4 interpolation nodes are: 

{(1,1), (2, 3), (4, 9), (8, 27)}. 

The vector with the interpolation data in the corresponding order is 


h= (0,3,43,345)^, 


and the coefficient matrix of the linear system is 


A 


All 1 \ 
13 2 6 
1 9 4 36 
^1 27 8 216y 


The solution of the linear system Ac = b is the vector 


c=(2,-3,-l,2)'^, 


and therefore the implicit equation of the curve is 


2 — 3y — X + 2xy = 0. 


In [13], where this approach is presented for the case of surface implicitization, 
it is indicated that in this way the problem reduces to interpolating a univari¬ 
ate polynomial. However, it must be observed that the coefficient matrix of 
the linear system is not a Vandermonde matrix associated with a univariate 
Lagrange interpolation problem; it is the transpose of such a matrix. 

The approach of [13] has its roots in computer algebra (see, for example, [23]), 
and so the solution of Vandermonde systems by Bjorck-Pereyra algorithms [1] 
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is not considered there. A recent extension of the Bjorck-Pereyra approach 
(closely related to bidiagonal factorizations of the inverse of the matrix) to 
Vandermonde-like matrices can be seen in Chapter 22 of [9]. 


4 A choice of nodes leading to the Kronecker product 


Another approach for computing the implicit equation of a plane rationally 
parametrized curve by means of resultants and interpolation is the one intro¬ 
duced in [15]. 

In that paper, the polynomial dehning the implicit equation of the curve, 
that is, Restiuiit) — xvi{t),U 2 {t) — yv 2 {t)), is computed by using a bivariate 
interpolation technique in which the nodes are arranged forming a tensor 
product grid. This choice of the interpolation nodes is specially appropriate 
for the interpolation space we are working with, because it reduces the solution 
of the interpolation problem to the solution of a linear system of order N whose 
coefficient matrix is the Kronecker product 

A = 14 (g) Vy, 

(where the Kronecker product B ^ D is dehned by blocks as {bkiD), with 
B = {bki)) with 14 being the Vandermonde matrix generated by the hrst 
component of the interpolation nodes and Vy being the Vandermonde matrix 
generated by the second component of the interpolation nodes (see [15] for the 
details). In addition, the algorithm included there reduces the solution of this 
linear system with Kronecker product structure to solving m-|-l Vandermonde 
linear systems with the same matrix Vy and n-|-l Vandermonde linear systems 
with the same matrix 14 - In this way, the solution of a bivariate interpolation 
problem is reduced to the solution of only univariate interpolation problems 
in the variables x and y. 

A Maple implementation of the complete algorithm is also included in the 
paper. As every linear system to be solved is a Vandermonde linear system, 
the algorithm uses the Bjorck-Pereyra algorithm [1, 8] for solving them, since 
it takes advantage of the special structure of the coefficient matrices 14 and 
V,. 

The specific choice of the interpolation nodes proposed in the paper is 

= (hi) : i = 0;...,n}, 

in the same order as the interpolation basis. 
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Now we apply this technique to our example. 

The interpolation nodes considered in this case are: 


{( 0 , 0 ),( 0 , 1 ),( 1 , 0 ),( 1 , 1 )}. 

The vector with the interpolation data in the corresponding order is: 

6 =( 2 ,- 1 , 1 , 0 )^. 

The coefficient matrix of the linear system is: 

^1 0 0 0 ^ 

110 0 
10 10 
^1 1 1 ly 

The solution of the linear system Ac = b is the vector 

c=(2,-3,-l,2)^, 

and therefore, the implicit equation of the curve is 

2 — 3y — X + 2xy = 0. 


^1 0^ 


\ oV 



V 


5 Some remarks on the computational complexity 


We start this section by presenting a bigger (but still small) example of curve 
implicitization in which we will show the different behaviour of the three 
different approaches we have described in the previous sections. 


Let 


P = {x{t),y{t)) 


/2f2 + 2t+l + 

V P+5 ’ P-3 ) 


be a rational parametrization of a curve C whose implicit equation is given 
by the polynomial 


F{x,y) = —53 + 4:2y — 74:y‘^ + 172x + 707xy + 121xy‘^ + 37xy^ — 652x‘^ — 1156x‘^y 
—490x^y^ — 34x‘^y^ + 626x^ + 396x^y + 432x^1/^ — 2x^y^. 
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When computing F{x^ y) by means of the approach described in Section 2, 
the implicitization problem is reduced to the solution of a homogeneous linear 
system of order 16 with a non-structured coefficient matrix A. As this matrix 
is too big for including it here, we just show one of its largest entries: 

^ _ 761421163154846949 

“ 149346877368718693 ■ 

After solving this linear system the vector with the coefficients of F{x,y) is 
obtained. 


The other two methods need the Sylvester matrix of the polynomials p{t) = 
2F + 2t + 1 — x{F + 5) and q{t) = F — SF + t — 1 — y{F — 3) for computing 
the implicit equation of C. This Sylvester matrix is: 

(-x 2 2 -bx +1 0 0 \ 


0 

—X 

2 

2 

—5x + 1 

0 

0 

0 

—X 

2 

2 

—5x -f-1 

1 

-y-3 

1 

3y-l 

0 

0 

0 

1 

-y-3 

1 

3y-l 

0 


0 

1 

-y-3 

1 

31/-1 y 


When using the approach described in Section 3, the implicitization problem 
is reduced to the solution of a non-homogeneous linear system Ac = 5 of 
order 16 where A is the transpose of a Vandermonde matrix. As we have done 
before, we only present the greatest entry of the matrix, which corresponds to 
the evaluation of the monomial x^y^ at the interpolation node (2^^, 3^^): 

Ai6,i6 = 103945637534048876111514866313854976. 

As for the vector b with the interpolation data, its largest component is 

bie = -207995995871362988895940143529893921. 

The solution of this linear system in which very large numbers are involved is 
the vector with the coefficients of F{x,y). 


Finally, when we compute the implicit equation of C by means of the method 
described in Section 4, the problem is reduced to the solution of the non- 
homogeneous linear system Ac = 5 of order 16 where A is the Kronecker 
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product of two Vandermonde matrices of order 4, 


^1 0 0 0^ 


^1 0 0 0^ 

1111 

0 

1111 

12 4 8 


12 4 8 

^1 3 9 27j 


^1 3 9 27j 


and the vector b containing the interpolation data is 

(- 53 , - 85 , - 265 , - 593 , 93 , 72 , 35 , - 12 , 2691 , 4277 , 8723 , 15561 , 11497 , 21242 , 44579 , 80014 )^. 

In this case only 8 linear systems with the Vandermonde matrix of order 4 
presented above have to be solved for obtaining the coefficients of F{x, y). 


Remark. Let us observe here that, as it can be read in [20], the implicit 
equation of a rational curve can also be computed by using techniques based 
on computing Grobner bases with the pure lexicographical ordering. However, 
although these techniques are very important from a theoretical point of view, 
they are not so effective in practice because they are very time and space 
consuming, even for problems of moderate degree. An example of this situation 
is the example introduced in this section. When we tried to compute the 
implicit equation of C by using the Maple command for computing Grobner 
bases gbasis no answer was obtained after a lot of minutes of computation, 
and in addition a lot of memory space (more that 70 megabytes without 
hnishing the computation) was required. 

In this sense, it can be read in [7] that the complexity theory for Grobner bases 
gives rise to the pessimistic view that these methods for polynomial ideals are 
not useful in practice, except for rather small cases. In [20] it is recognized that 
the use of Grobner bases for surface implicitization is not very computationally 
efficient, and this fact is also observed in a different application of Grobner 
bases in [3]. As illustrated with our example in this section, the problems with 
Grobner bases appear also in the case of curve implicitization and with small 
degrees. 


As for the approaches based on interpolation, it must be observed that the 
numbers involved in the hrst and in the second approach are much larger than 
the numbers involved in the third one (see the example at the beginning of 
this section), which can make the computations slower. 

Now we briefly analyze the computational complexity of the three methods 
for curve implicitization presented in this paper, beginning with the stage 
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corresponding to the solution of the linear systems. We will do it in terms of 
arithmetic operations, and for the sake of simplicity we assume m = n. 

- The computational complexity of the first approach is the computational 
complexity of solving a non-structured linear system of order N: 0{N^). 

- The computational complexity of the second approach is the computational 
complexity of solving a linear system whose coefficient matrix is the trans¬ 
pose of a Vandermonde matrix of order N: 0{N‘^) if the Bjorck-Pereyra 
algorithm is used [1, 8]. 

- The computational complexity of the third approach is the computational 

complexity of solving 2{n+ 1) Vandermonde linear systems of order (n-|- 1): 
0{{n -|- 1)^) = if the Bjorck-Pereyra algorithm is used for solving 

each Vandermonde linear system [1, 8]. 

Let us point out that in the first approach, in addition to the higher compu¬ 
tational cost of solving the linear system, hrst it is necessary to evaluate the 
parametric equations to obtain the interpolation nodes, and then the coeffi¬ 
cient matrix of the linear system (which is not a structured matrix) must be 
constructed and stored. 

As for the computation of the interpolation data, its complexity is the same 
both in the second and in the third approach, and it is the complexity of 
computing (n -|- 1)^ determinants of order 0{n): 0{n^) = Let us 

point out here that, although this complexity is greater than the complexity 
of solving the linear systems, this stage can easily be parallelized because each 
datum can be computed separately. 

The design of parallel algorithms for computer algebra problems such as re¬ 
sultant computation and elimination of variables has been recently considered 
in [2] and [10]. In this sense, let us observe that the approach introduced in 
Section 4 has a high degree of intrinsic parallelism present not only in the 
computation of the interpolation data but also in the computation of the co¬ 
efficients of the interpolating polynomial, since each one of the linear systems 
with matrix Vy can be solved independently and the same happens to each 
one of the linear systems with matrix 14 (see [15]). 

Although it is true that this measure of computational complexity -which is 
the standard one in numerical computations- is not a complete measure of 
the computational complexity in exact arithmetic, where the computational 
cost also depends on the size and structure of the numbers, it nevertheless can 
serve to have useful estimates of the computational complexity. In this sense, 
we can read in [19] how the reduction of the order of the resultant matrix from 
2n to n may lead to faster computations. 

In addition, it must be pointed out that the two implicitization methods based 
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on interpolation and resultants we have described can also take advantage of 
this reduction in the order of the resultant matrix because they can be used 
with any other resultant matrix, for example, the resultant matrix obtained 
when using the method of moving curves [19]. 


Finally, it must be noticed that in the situation described in the approach of 
Section 4, that is, when the nodes are arranged forming a Cartesian product 
grid, explicit formulas (involving the Lagrange basis) exist for computing the 
interpolation polynomial [18]. 

Nevertheless, the existence of an explicit formula does not imply there is no 
computational cost in applying it. In this sense we can recall the classical 
paper [6], where it is stressed that the fact that we have Cramer’s rule does 
not make the practical solution of linear systems a trivial and dull task. 

In our case, it must be taken into account that the computation of the inter¬ 
polation polynomial in the monomial basis by using an explicit formula has 
a computational cost (see [7] for the cost in the univariate case). It must be 
observed that our approach described in Section 4 has a complexity of the 
same order as the algorithm presented in a more general setting in [16]. 
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